#include<TELFile.h>

#include<TCanvas.h>
#include<TH1F.h>
#include<TH2F.h>
#include<TF1.h>

TELFile::TELFile(Int_t RN)
{

  RunNumber = RN;

  glob = new Global(RunNumber);

  Int_t StringLength = 200;

  SiFilePath = new Char_t[StringLength];
  SiFilePath = glob->TELRawPath;

  SiFile = fopen(SiFilePath,"r");

  if(!SiFile){
    cout << "!!! ERROR !!!" << endl;
    cout << "File " << SiFilePath << " does not exist" << endl;
    return;
  }

  ADCXSize = 0;
  ADCYSize = 0;

  ShiftX = new Int_t[glob->TELNumberOfPlanes];
  ShiftY = new Int_t[glob->TELNumberOfPlanes];

  for(Int_t i = 0; i < glob->TELNumberOfPlanes; i++){
    ShiftX[i] = 0;
    ShiftY[i] = 0;
  }

  for(Int_t i = 0; i < glob->TELNumberOfPlanes; i++){
    ADCXSize += glob->TELNumberOfStripsX[i];
    ADCYSize += glob->TELNumberOfStripsY[i];
    if(i > 0){
      ShiftX[i] = ShiftX[i-1] + glob->TELNumberOfStripsX[i-1];
      ShiftY[i] = ShiftY[i-1] + glob->TELNumberOfStripsY[i-1];
    }
 }

  PlaneXADC = new Int_t[ADCXSize];
  PlaneYADC = new Int_t[ADCYSize];

  TELRoot = new Char_t[StringLength];
  sprintf(TELRoot,"ntuple/tele/TELData%u.root",RunNumber);
  
  TELData = new TFile(TELRoot);

  if(TELData->IsZombie()){

    TELData = new TFile(TELRoot,"RECREATE");
    
    TELTree = new TTree("TELData","Telescope Data");
    TELTree->Branch("Trigger",&Trigger,"Trigger[2]/D");
    
    sprintf(TELRoot,"PlaneXADC[%u]/I",ADCXSize);
    TELTree->Branch("PlaneX",&PlaneXADC[0],TELRoot);
    sprintf(TELRoot,"PlaneYADC[%u]/I",ADCYSize);
    TELTree->Branch("PlaneY",&PlaneYADC[0],TELRoot);
    
    //Run Header
    Char_t Header[80];
    for(Int_t i = 0; i < 5; i++){
      fgets(Header,80,SiFile); printf("%s",Header);
    }
    Int_t eventcounter = 0;
    while(ReadEvents()) eventcounter++;
    cout << "Number of Events in file: " << eventcounter << endl;
    TELData->cd();
    TELTree->Print();
    TELTree->Write();
    TELData->Close();

  }

  APara = new Double_t[(glob->TELNumberOfPlanes - 1)*6];

}//TELFile

Int_t TELFile::ReadEvents(){
  Int_t TelPointer = 0;
  Int_t tmp;
  Int_t ADC[2];

  for(Int_t i = 0; i < ADCXSize; i++) PlaneXADC[i] = 0;
  for(Int_t i = 0; i < ADCYSize; i++) PlaneYADC[i] = 0;

  for(Int_t i = 0; i < EventLength; i++){
    
    if(feof(SiFile) || ferror(SiFile)) return 0;
    
    fread(&tmp,4,1,SiFile);
    tmp = swap4(tmp);

    if(i == 0){
      if(tmp == 0xDDCCBBAA) break;
      cout << " -- Event -- " << tmp << endl;
    }

    if(i == 1) Trigger[0] = tmp;
    if(i == 2){
      Trigger[1] = tmp;
      cout << "Trigger time: " << Trigger[0] << " s " << Trigger[1] << " us" << endl; 
    }

    if(i >= 4){
     
      if(tmp == 0xAABBCCDD) break;
 
      //end of first telescope
      if(tmp == 0xAABBCC00){
	TelPointer = 1;
	continue;
      }

      //end of second telescope
      if(tmp == 0xAABBCC01){
	TelPointer = 2;
	continue;
      }

      //end of second telescope
      if(tmp == 0xAABBCC03) continue;

      ADC[0] = (tmp >> 12)&0xfff;
      ADC[1] = (tmp & 0xfff);

      if(ADC[0] < glob->TELNumberOfStripsX[TelPointer]){
	PlaneXADC[ADC[0] + ShiftX[TelPointer]] = ADC[1];
      }
      
      if(ADC[0] >= glob->TELNumberOfStripsX[TelPointer]){
	PlaneYADC[ADC[0] - glob->TELNumberOfStripsX[TelPointer] + ShiftY[TelPointer]] = ADC[1];
      }

    }//if(i >= 4)

  }//for

  if(!feof(SiFile)) TELTree->Fill();
  if(feof(SiFile) || ferror(SiFile)) return 0;

  return 1;

}

Int_t TELFile::swap4(Int_t w4){

  unsigned int rw4;
  unsigned char b[4];

  for(Int_t i = 0; i < 4; i++) b[i] = (w4>>i*8)&0xff; 

  rw4 = 0; 

  for(Int_t i = 0; i < 4; i++) rw4+=b[i]<<(3-i)*8;
 
  return (Int_t)rw4;

}

void TELFile::SimplePedNoi(){

  PedX = new Float_t[ADCXSize];
  PedY = new Float_t[ADCYSize];
  NoiX = new Float_t[ADCXSize];
  NoiY = new Float_t[ADCYSize];

  Float_t* PedX2 = new Float_t[ADCXSize];
  Float_t* NormX = new Float_t[ADCXSize];
  Float_t* PedY2 = new Float_t[ADCYSize];
  Float_t* NormY = new Float_t[ADCYSize];

  Float_t HistMaxX[glob->TELNumberOfPlanes];
  Float_t HistMaxY[glob->TELNumberOfPlanes];

  for(Int_t i = 0; i < glob->TELNumberOfPlanes; i++){
    HistMaxX[i] = 0.0;
    HistMaxY[i] = 0.0;
  }

  for(Int_t i_ped = 0; i_ped < ADCXSize; i_ped++){
    PedX[i_ped] = 0.0;
    NoiX[i_ped] = 0.0;
    PedX2[i_ped] = 0.0;
    NormX[i_ped] = 0.0;
  }

  for(Int_t i_ped = 0; i_ped < ADCYSize; i_ped++){
    PedY[i_ped] = 0.0;
    NoiY[i_ped] = 0.0;
    PedY2[i_ped] = 0.0;
    NormY[i_ped] = 0.0;
  }

  //Histograms for pedestal
  TH1F* h11[2][glob->TELNumberOfPlanes];
  TCanvas* canv11 = new TCanvas("c1","",1200,800);
  canv11->SetFillColor(0);
  canv11->Divide(glob->TELNumberOfPlanes,2);

  //Histograms for noise
  TH1F* h12[2][glob->TELNumberOfPlanes];
  TCanvas* canv12 = new TCanvas("c2","",1200,800);
  canv12->SetFillColor(0);
  canv12->Divide(glob->TELNumberOfPlanes,2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    sprintf(TELRoot,"Ped_X_%u",n+1);
    h11[0][n] = new TH1F(TELRoot,"",glob->TELNumberOfStripsX[n],0,glob->TELNumberOfStripsX[n]); 
    h11[0][n]->GetXaxis()->SetTitle("Chanel number");
    h11[0][n]->GetYaxis()->SetTitle("Pedestal [ADC]");
    sprintf(TELRoot,"Ped_X_%u",n+1);
    h11[0][n]->SetTitle(TELRoot);
    sprintf(TELRoot,"Ped_Y_%u",n+1);
    h11[1][n] = new TH1F(TELRoot,"",glob->TELNumberOfStripsY[n],0,glob->TELNumberOfStripsY[n]);
    h11[1][n]->GetXaxis()->SetTitle("Chanel number");
    h11[1][n]->GetYaxis()->SetTitle("Pedestal [ADC]");
    sprintf(TELRoot,"Ped_Y_%u",n+1);
    h11[1][n]->SetTitle(TELRoot);

    sprintf(TELRoot,"Noi_X_%u",n+1);
    h12[0][n] = new TH1F(TELRoot,"",glob->TELNumberOfStripsX[n],0,glob->TELNumberOfStripsX[n]); 
    h12[0][n]->GetXaxis()->SetTitle("Chanel number");
    h12[0][n]->GetYaxis()->SetTitle("Noise [ADC]");
    sprintf(TELRoot,"Noi_X_%u",n+1);
    h12[0][n]->SetTitle(TELRoot);
    sprintf(TELRoot,"Noi_Y_%u",n+1);
    h12[1][n] = new TH1F(TELRoot,"",glob->TELNumberOfStripsY[n],0,glob->TELNumberOfStripsY[n]);
    h12[1][n]->GetXaxis()->SetTitle("Chanel number");
    h12[1][n]->GetYaxis()->SetTitle("Noise [ADC]");
    sprintf(TELRoot,"Noi_Y_%u",n+1);
    h12[1][n]->SetTitle(TELRoot);

  }

  sprintf(TELRoot,"ntuple/tele/TELData%u.root",RunNumber);

  TELData = new TFile(TELRoot);
  TTree* ntuple = (TTree*)TELData->Get("TELData");

  ntuple->SetBranchAddress("PlaneX",&PlaneXADC[0]);
  ntuple->SetBranchAddress("PlaneY",&PlaneYADC[0]);  

  for(Int_t i = 0; i < ntuple->GetEntries(); i++){

    ntuple->GetEntry(i);

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

      for(Int_t x = 0; x < glob->TELNumberOfStripsX[n]; x++){
	
	PedX[x + ShiftX[n]] += (Float_t)PlaneXADC[x + ShiftX[n]];
	PedX2[x + ShiftX[n]] += TMath::Power((Float_t)PlaneXADC[x + ShiftX[n]],2);
	NormX[x + ShiftX[n]] += 1.0;
	
      }//x

      for(Int_t y = 0; y < glob->TELNumberOfStripsY[n]; y++){
	
	PedY[y + ShiftY[n]] += (Float_t)PlaneYADC[y + ShiftY[n]];
	PedY2[y + ShiftY[n]] += TMath::Power((Float_t)PlaneYADC[y + ShiftY[n]],2);
        NormY[y + ShiftY[n]] += 1.0;
	
      }//y
      
    }//n

  }//i

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    for(Int_t x = 0; x < glob->TELNumberOfStripsX[n]; x++){
      if(NormX[x + ShiftX[n]] <= 1.0){
	PedX[x + ShiftX[n]] = 0.0;
	NoiX[x + ShiftX[n]] = 0.0;
	h11[0][n]->Fill(x,PedX[x + ShiftX[n]]);
	h12[0][n]->Fill(x,NoiX[x + ShiftX[n]]);
	continue;
      }
      PedX[x + ShiftX[n]] = PedX[x + ShiftX[n]]/NormX[x + ShiftX[n]];
      PedX2[x + ShiftX[n]] = PedX2[x + ShiftX[n]]/NormX[x + ShiftX[n]];
      NoiX[x + ShiftX[n]] = TMath::Sqrt((NormX[x + ShiftX[n]]/(NormX[x + ShiftX[n]] - 1))*(PedX2[x + ShiftX[n]] - TMath::Power(PedX[x + ShiftX[n]],2)));
      h11[0][n]->Fill(x,PedX[x + ShiftX[n]]);
      h12[0][n]->Fill(x,NoiX[x + ShiftX[n]]);
    }//x

    for(Int_t y = 0; y < glob->TELNumberOfStripsY[n]; y++){
      if(NormY[y + ShiftY[n]] <= 1.0){
	PedY[y + ShiftY[n]] = 0.0;
	NoiY[y + ShiftY[n]] = 0.0;
	h11[1][n]->Fill(y,PedY[y + ShiftY[n]]);
	h12[1][n]->Fill(y,NoiY[y + ShiftY[n]]);
	continue;
      }
      PedY[y + ShiftY[n]] = PedY[y + ShiftY[n]]/NormY[y + ShiftY[n]];
      PedY2[y + ShiftY[n]] = PedY2[y + ShiftY[n]]/NormY[y + ShiftY[n]];
      NoiY[y + ShiftY[n]] = TMath::Sqrt((NormY[y + ShiftY[n]]/(NormY[y + ShiftY[n]] - 1))*(PedY2[y + ShiftY[n]] - TMath::Power(PedY[y + ShiftY[n]],2)));
      h11[1][n]->Fill(y,PedY[y + ShiftY[n]]);
      h12[1][n]->Fill(y,NoiY[y + ShiftY[n]]);
    }

  }//n

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    canv11->cd(n+1);
    h11[0][n]->Draw();
    canv11->cd(n + 1  + glob->TELNumberOfPlanes);
    h11[1][n]->Draw();

    canv12->cd(n+1);
    h12[0][n]->Draw();
    canv12->cd(n + 1 + glob->TELNumberOfPlanes);
    h12[1][n]->Draw();

  }

  //Histograms for signals
  TH1F* h13[2][glob->TELNumberOfPlanes];
  TCanvas* canv13 = new TCanvas("c3","",1200,800);
  canv13->SetFillColor(0);
  canv13->Divide(glob->TELNumberOfPlanes,2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
    sprintf(TELRoot,"Signal_X_%u",n+1);
    h13[0][n] = new TH1F(TELRoot,"",100,0,HistMaxX[n]); 
    //h13[0][n]->GetXaxis()->SetTitle("Cluster charge [ADC]");
    h13[0][n]->GetXaxis()->SetTitle("Pedestals [ADC]"); 
    h13[0][n]->GetYaxis()->SetTitle("Number of events");
    sprintf(TELRoot,"Telescope %u (plane X)",n+1);
    h13[0][n]->SetTitle(TELRoot);
    sprintf(TELRoot,"Signal_Y_%u",n+1);
    h13[1][n] = new TH1F(TELRoot,"",100,0,HistMaxY[n]);
    //    h13[1][n]->GetXaxis()->SetTitle("Cluster charge [ADC]");
    h13[1][n]->GetXaxis()->SetTitle("Pedestal [ADC]");
    h13[1][n]->GetYaxis()->SetTitle("Number of events");
    sprintf(TELRoot,"Telescope %u (plane Y)",n+1);
    h13[1][n]->SetTitle(TELRoot);
  }

  for(Int_t i = 0; i < ntuple->GetEntries(); i++){

    ntuple->GetEntry(i);

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

      for(Int_t x = 0; x < glob->TELNumberOfStripsX[n]; x++){

	if(PlaneXADC[x + ShiftX[n]] > 0) h13[0][n]->Fill((Float_t)PlaneXADC[x + ShiftX[n]]);

      }//x

      for(Int_t y = 0; y < glob->TELNumberOfStripsY[n]; y++){

	if(PlaneYADC[y + ShiftY[n]] > 0) h13[1][n]->Fill((Float_t)PlaneYADC[y + ShiftY[n]]);

      }//y

    }//n

  }//i

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    canv13->cd(n+1);
    h13[0][n]->Draw();
    canv13->cd(n + 1  + glob->TELNumberOfPlanes);
    h13[1][n]->Draw();

  }

  delete [] PedX;
  delete [] PedY;
  delete [] NoiX;
  delete [] NoiY;

  delete [] PedX2;
  delete [] PedY2;
  delete [] NormX;
  delete [] NormY;

}//SimplePedNoi

void TELFile::PedNoi(){

  PedX = new Float_t[ADCXSize];
  PedY = new Float_t[ADCYSize];
  NoiX = new Float_t[ADCXSize];
  NoiY = new Float_t[ADCYSize];

  Float_t* PedX2 = new Float_t[ADCXSize];
  Float_t* NormX = new Float_t[ADCXSize];
  Float_t* PedY2 = new Float_t[ADCYSize];
  Float_t* NormY = new Float_t[ADCYSize];

  Float_t HistMaxX[glob->TELNumberOfPlanes];
  Float_t HistMaxY[glob->TELNumberOfPlanes];

  for(Int_t i = 0; i < glob->TELNumberOfPlanes; i++){
    HistMaxX[i] = 0.0;
    HistMaxY[i] = 0.0;
  }

  for(Int_t i_ped = 0; i_ped < ADCXSize; i_ped++){
    PedX[i_ped] = 0.0;
    NoiX[i_ped] = 0.0;
    PedX2[i_ped] = 0.0;
    NormX[i_ped] = 0.0;
  }

  for(Int_t i_ped = 0; i_ped < ADCYSize; i_ped++){
    PedY[i_ped] = 0.0;
    NoiY[i_ped] = 0.0;
    PedY2[i_ped] = 0.0;
    NormY[i_ped] = 0.0;
  }

  //Histograms for pedestal
  TH1F* h11[2][glob->TELNumberOfPlanes];
  TCanvas* canv11 = new TCanvas("c1","",1200,800);
  canv11->SetFillColor(0);
  canv11->Divide(glob->TELNumberOfPlanes,2);

  //Histograms for noise
  TH1F* h12[2][glob->TELNumberOfPlanes];
  TCanvas* canv12 = new TCanvas("c2","",1200,800);
  canv12->SetFillColor(0);
  canv12->Divide(glob->TELNumberOfPlanes,2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    sprintf(TELRoot,"Ped_X_%u",n+1);
    h11[0][n] = new TH1F(TELRoot,"",glob->TELNumberOfStripsX[n],0,glob->TELNumberOfStripsX[n]); 
    h11[0][n]->GetXaxis()->SetTitle("Chanel number");
    h11[0][n]->GetYaxis()->SetTitle("Pedestal [ADC]");
    sprintf(TELRoot,"Ped_X_%u",n+1);
    h11[0][n]->SetTitle(TELRoot);
    sprintf(TELRoot,"Ped_Y_%u",n+1);
    h11[1][n] = new TH1F(TELRoot,"",glob->TELNumberOfStripsY[n],0,glob->TELNumberOfStripsY[n]);
    h11[1][n]->GetXaxis()->SetTitle("Chanel number");
    h11[1][n]->GetYaxis()->SetTitle("Pedestal [ADC]");
    sprintf(TELRoot,"Ped_Y_%u",n+1);
    h11[1][n]->SetTitle(TELRoot);

    sprintf(TELRoot,"Noi_X_%u",n+1);
    h12[0][n] = new TH1F(TELRoot,"",glob->TELNumberOfStripsX[n],0,glob->TELNumberOfStripsX[n]); 
    h12[0][n]->GetXaxis()->SetTitle("Chanel number");
    h12[0][n]->GetYaxis()->SetTitle("Noise [ADC]");
    sprintf(TELRoot,"Noi_X_%u",n+1);
    h12[0][n]->SetTitle(TELRoot);
    sprintf(TELRoot,"Noi_Y_%u",n+1);
    h12[1][n] = new TH1F(TELRoot,"",glob->TELNumberOfStripsY[n],0,glob->TELNumberOfStripsY[n]);
    h12[1][n]->GetXaxis()->SetTitle("Chanel number");
    h12[1][n]->GetYaxis()->SetTitle("Noise [ADC]");
    sprintf(TELRoot,"Noi_Y_%u",n+1);
    h12[1][n]->SetTitle(TELRoot);

  }

  Int_t HRejectMax = 200;
  TH1F* HReject;

  sprintf(TELRoot,"ntuple/tele/TELData%u.root",RunNumber);

  TELData = new TFile(TELRoot);
  TTree* ntuple = (TTree*)TELData->Get("TELData");

  ntuple->SetBranchAddress("PlaneX",&PlaneXADC[0]);
  ntuple->SetBranchAddress("PlaneY",&PlaneYADC[0]);  

  for(Int_t i = 0; i < ntuple->GetEntries(); i++){

    ntuple->GetEntry(i);

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

      HReject = new TH1F("HReject","",HRejectMax,0,HRejectMax);

      for(Int_t x = 0; x < glob->TELNumberOfStripsX[n]; x++) if(PlaneXADC[x + ShiftX[n]] < HRejectMax) HReject->Fill((Float_t)PlaneXADC[x + ShiftX[n]]); 

      for(Int_t x = 0; x < glob->TELNumberOfStripsX[n]; x++){

	if(HistMaxX[n] < (Float_t)PlaneXADC[x + ShiftX[n]]) HistMaxX[n] = (Float_t)PlaneXADC[x + ShiftX[n]];  

	if((Float_t)PlaneXADC[x + ShiftX[n]] > HReject->GetMean() + 3.0*HReject->GetRMS()) continue;

	PedX[x + ShiftX[n]] += (Float_t)PlaneXADC[x + ShiftX[n]];
	PedX2[x + ShiftX[n]] += TMath::Power((Float_t)PlaneXADC[x + ShiftX[n]],2);
	NormX[x + ShiftX[n]] += 1.0;

      }//x

      delete HReject;
      HReject = new TH1F("HReject","",HRejectMax,0,HRejectMax);

      for(Int_t y = 0; y < glob->TELNumberOfStripsY[n]; y++) if(PlaneYADC[y + ShiftY[n]] < HRejectMax) HReject->Fill((Float_t)PlaneYADC[y + ShiftY[n]]);

      for(Int_t y = 0; y < glob->TELNumberOfStripsY[n]; y++){

        if(HistMaxY[n] < (Float_t)PlaneYADC[y + ShiftY[n]]) HistMaxY[n] = (Float_t)PlaneYADC[y + ShiftY[n]];

        if((Float_t)PlaneYADC[y + ShiftY[n]] > HReject->GetMean() + 3.0*HReject->GetRMS()) continue;

	PedY[y + ShiftY[n]] += (Float_t)PlaneYADC[y + ShiftY[n]];
	PedY2[y + ShiftY[n]] += TMath::Power((Float_t)PlaneYADC[y + ShiftY[n]],2);
        NormY[y + ShiftY[n]] += 1.0;

      }//y

      delete HReject;

    }//n

  }//i

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    for(Int_t x = 0; x < glob->TELNumberOfStripsX[n]; x++){
      if(NormX[x + ShiftX[n]] <= 1.0){
	PedX[x + ShiftX[n]] = 0.0;
	NoiX[x + ShiftX[n]] = 0.0;
	h11[0][n]->Fill(x,PedX[x + ShiftX[n]]);
	h12[0][n]->Fill(x,NoiX[x + ShiftX[n]]);
	continue;
      }
      PedX[x + ShiftX[n]] = PedX[x + ShiftX[n]]/NormX[x + ShiftX[n]];
      PedX2[x + ShiftX[n]] = PedX2[x + ShiftX[n]]/NormX[x + ShiftX[n]];
      NoiX[x + ShiftX[n]] = TMath::Sqrt((NormX[x + ShiftX[n]]/(NormX[x + ShiftX[n]] - 1))*(PedX2[x + ShiftX[n]] - TMath::Power(PedX[x + ShiftX[n]],2)));
      h11[0][n]->Fill(x,PedX[x + ShiftX[n]]);
      h12[0][n]->Fill(x,NoiX[x + ShiftX[n]]);
    }//x

    for(Int_t y = 0; y < glob->TELNumberOfStripsY[n]; y++){
      if(NormY[y + ShiftY[n]] <= 1.0){
	PedY[y + ShiftY[n]] = 0.0;
	NoiY[y + ShiftY[n]] = 0.0;
	h11[1][n]->Fill(y,PedY[y + ShiftY[n]]);
	h12[1][n]->Fill(y,NoiY[y + ShiftY[n]]);
	continue;
      }
      PedY[y + ShiftY[n]] = PedY[y + ShiftY[n]]/NormY[y + ShiftY[n]];
      PedY2[y + ShiftY[n]] = PedY2[y + ShiftY[n]]/NormY[y + ShiftY[n]];
      NoiY[y + ShiftY[n]] = TMath::Sqrt((NormY[y + ShiftY[n]]/(NormY[y + ShiftY[n]] - 1))*(PedY2[y + ShiftY[n]] - TMath::Power(PedY[y + ShiftY[n]],2)));
      h11[1][n]->Fill(y,PedY[y + ShiftY[n]]);
      h12[1][n]->Fill(y,NoiY[y + ShiftY[n]]);
    }

  }//n

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    canv11->cd(n+1);
    h11[0][n]->Draw();
    canv11->cd(n + 1  + glob->TELNumberOfPlanes);
    h11[1][n]->Draw();

    canv12->cd(n+1);
    h12[0][n]->Draw();
    canv12->cd(n + 1 + glob->TELNumberOfPlanes);
    h12[1][n]->Draw();

  }

  //Histograms for signals
  TH1F* h13[2][glob->TELNumberOfPlanes];
  TCanvas* canv13 = new TCanvas("c3","",1200,800);
  canv13->SetFillColor(0);
  canv13->Divide(glob->TELNumberOfPlanes,2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
    sprintf(TELRoot,"Signal_X_%u",n+1);
    h13[0][n] = new TH1F(TELRoot,"",100,0,HistMaxX[n]); 
    h13[0][n]->GetXaxis()->SetTitle("Cluster charge [ADC]");
    h13[0][n]->GetYaxis()->SetTitle("Number of events");
    sprintf(TELRoot,"Telescope %u (plane X)",n+1);
    h13[0][n]->SetTitle(TELRoot);
    sprintf(TELRoot,"Signal_Y_%u",n+1);
    h13[1][n] = new TH1F(TELRoot,"",100,0,HistMaxY[n]);
    h13[1][n]->GetXaxis()->SetTitle("Cluster charge [ADC]");
    h13[1][n]->GetYaxis()->SetTitle("Number of events");
    sprintf(TELRoot,"Telescope %u (plane Y)",n+1);
    h13[1][n]->SetTitle(TELRoot);
  }

  for(Int_t i = 0; i < ntuple->GetEntries(); i++){

    ntuple->GetEntry(i);

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

      for(Int_t x = 0; x < glob->TELNumberOfStripsX[n]; x++){

	if(PlaneXADC[x + ShiftX[n]] > 0) h13[0][n]->Fill((Float_t)PlaneXADC[x + ShiftX[n]]);

      }//x

      for(Int_t y = 0; y < glob->TELNumberOfStripsY[n]; y++){

	if(PlaneYADC[y + ShiftY[n]] > 0) h13[1][n]->Fill((Float_t)PlaneYADC[y + ShiftY[n]]);

      }//y

    }//n

  }//i

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    canv13->cd(n+1);
    h13[0][n]->Draw();
    canv13->cd(n + 1  + glob->TELNumberOfPlanes);
    h13[1][n]->Draw();

  }

  delete [] PedX;
  delete [] PedY;
  delete [] NoiX;
  delete [] NoiY;

  delete [] PedX2;
  delete [] PedY2;
  delete [] NormX;
  delete [] NormY;

}//PedNoi

void TELFile::FindClusters(Int_t ClustSizeX, Int_t ClustSizeY){

  Float_t DigX[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Float_t DigY[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Float_t CogX[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Float_t CogY[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Float_t TrigTime[2];

  vector<Double_t> VecDigX[glob->TELNumberOfPlanes - 1][glob->TELNumberOfPlanes];
  vector<Double_t> VecDigY[glob->TELNumberOfPlanes - 1][glob->TELNumberOfPlanes];
  vector<Double_t> VecCogX[glob->TELNumberOfPlanes - 1][glob->TELNumberOfPlanes];
  vector<Double_t> VecCogY[glob->TELNumberOfPlanes - 1][glob->TELNumberOfPlanes];
  Bool_t VecFill;

  Bool_t* SeedX = new Bool_t[ADCXSize];
  Bool_t* SeedY = new Bool_t[ADCYSize];

  Float_t* HitBufferX = new Float_t[ADCXSize];
  Float_t* HitBufferY = new Float_t[ADCYSize];
  Float_t* CogBufferX = new Float_t[ADCXSize];
  Float_t* CogBufferY = new Float_t[ADCYSize];

  sprintf(TELRoot,"ntuple/tele/TELData%u.root",RunNumber);

  TELData = new TFile(TELRoot);
  TTree* ntuple = (TTree*)TELData->Get("TELData");

  ntuple->SetBranchAddress("PlaneX",&PlaneXADC[0]);
  ntuple->SetBranchAddress("PlaneY",&PlaneYADC[0]);  

  Float_t ClustCharge;

  //Histograms for signals
  TH1F* h14[2][glob->TELNumberOfPlanes];
  TCanvas* canv14 = new TCanvas("c4","",1200,800);
  canv14->SetFillColor(0);
  canv14->Divide(glob->TELNumberOfPlanes,2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
    sprintf(TELRoot,"Signal_X_%u",n+1);
    h14[0][n] = new TH1F(TELRoot,"",100,0,3000); 
    h14[0][n]->GetXaxis()->SetTitle("Cluster charge [ADC]");
    h14[0][n]->GetYaxis()->SetTitle("Number of events");
    sprintf(TELRoot,"Telescope %u (plane X)",n+1);
    h14[0][n]->SetTitle(TELRoot);
    h14[0][n]->SetLineWidth(3);
    h14[0][n]->SetLineColor(4);
    h14[0][n]->GetXaxis()->SetTitleSize(0.05);
    h14[0][n]->GetXaxis()->SetTitleOffset(0.9);

    sprintf(TELRoot,"Signal_Y_%u",n+1);
    h14[1][n] = new TH1F(TELRoot,"",100,0,3000);
    h14[1][n]->GetXaxis()->SetTitle("Cluster charge [ADC]");
    h14[1][n]->GetYaxis()->SetTitle("Number of events");
    sprintf(TELRoot,"Telescope %u (plane Y)",n+1);
    h14[1][n]->SetTitle(TELRoot);
    h14[1][n]->SetLineWidth(3);
    h14[1][n]->SetLineColor(4);
    h14[1][n]->GetXaxis()->SetTitleSize(0.05);
    h14[1][n]->GetXaxis()->SetTitleOffset(0.9);

  }

  Int_t NumberOfHitsX;
  Int_t NumberOfHitsY;

  TH1F* h15[2][glob->TELNumberOfPlanes];
  TCanvas* canv15 = new TCanvas("c5","",1200,800);
  canv15->SetFillColor(0);
  canv15->Divide(glob->TELNumberOfPlanes,2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
    sprintf(TELRoot,"NumX_%u",n+1);
    h15[0][n] = new TH1F(TELRoot,"",glob->TELMaxHit,0,glob->TELMaxHit); 
    h15[0][n]->GetXaxis()->SetTitle("Number of hits");
    h15[0][n]->GetYaxis()->SetTitle("Number of events");
    sprintf(TELRoot,"X_%u",n+1);
    h15[0][n]->SetTitle(TELRoot);
    sprintf(TELRoot,"NumY_%u",n+1);
    h15[1][n] = new TH1F(TELRoot,"",glob->TELMaxHit,0,glob->TELMaxHit);
    h15[1][n]->GetXaxis()->SetTitle("Number of hits");
    h15[1][n]->GetYaxis()->SetTitle("Number of events");
    sprintf(TELRoot,"Y_%u",n+1);
    h15[1][n]->SetTitle(TELRoot);
  }

  TH2F* h16[glob->TELNumberOfPlanes];
  TCanvas* canv16 = new TCanvas("c6","",1200,400);
  canv16->SetFillColor(0);
  canv16->Divide(glob->TELNumberOfPlanes);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    sprintf(TELRoot,"HitsX_%u",n+1);
    h16[n] = new TH2F(TELRoot,"",glob->TELNumberOfStripsX[n],0,glob->TELNumberOfStripsX[n] * glob->TELStripSizeX[n],glob->TELNumberOfStripsY[n],0,glob->TELNumberOfStripsY[n] * glob->TELStripSizeY[n]); 
    h16[n]->GetXaxis()->SetTitle("X [Chanel]");
    h16[n]->GetYaxis()->SetTitle("Y [Chanel]");
    sprintf(TELRoot,"X_%u",n+1);
    h16[n]->SetTitle(TELRoot);

  }

  //Histograms for noise
  TH2F* h17[2][glob->TELNumberOfPlanes - 1];
  TCanvas* canv17 = new TCanvas("c7","",600,600);
  canv17->SetFillColor(0);
  canv17->Divide(glob->TELNumberOfPlanes - 1,2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){

    sprintf(TELRoot,"CorrX_%u",n+1);
    h17[0][n] = new TH2F(TELRoot,"",glob->TELNumberOfStripsX[n],0,glob->TELNumberOfStripsX[n] * glob->TELStripSizeX[n],glob->TELNumberOfStripsX[n+1],0,glob->TELNumberOfStripsX[n+1] * glob->TELStripSizeX[n+1]); 
    sprintf(TELRoot,"X telescope %d",n+1);
    h17[0][n]->GetXaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"X telescope %d",n+2);
    h17[0][n]->GetYaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Correlation X_%u",n+1);
    h17[0][n]->SetTitle(TELRoot);

    sprintf(TELRoot,"CorrY_%u",n+1);
    h17[1][n] = new TH2F(TELRoot,"",glob->TELNumberOfStripsY[n],0,glob->TELNumberOfStripsY[n] * glob->TELStripSizeY[n],glob->TELNumberOfStripsY[n+1],0,glob->TELNumberOfStripsY[n+1] * glob->TELStripSizeY[n+1]); 
    sprintf(TELRoot,"Y telescope %d",n+1);
    h17[1][n]->GetXaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Y telescope %d",n+2);
    h17[1][n]->GetYaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Correlation Y_%u",n+1);
    h17[1][n]->SetTitle(TELRoot);

  }

  //Histograms for Cog
  TH1F* h18[2][glob->TELNumberOfPlanes];
  TCanvas* canv18 = new TCanvas("c8","",1200,800);
  canv18->SetFillColor(0);
  canv18->Divide(glob->TELNumberOfPlanes,2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
    sprintf(TELRoot,"Cog_X_%u",n+1);
    h18[0][n] = new TH1F(TELRoot,"",100,-2.0*glob->TELStripSizeX[n],2.0*glob->TELStripSizeX[n]); 
    h18[0][n]->GetXaxis()->SetTitle("X_{Cog}-X_{Dig}");
    h18[0][n]->GetYaxis()->SetTitle("Number of events");
    sprintf(TELRoot,"Cog_X_%u",n+1);
    h18[0][n]->SetTitle(TELRoot);
    sprintf(TELRoot,"Cog_Y_%u",n+1);
    h18[1][n] = new TH1F(TELRoot,"",100,-2.0*glob->TELStripSizeY[n],2.0*glob->TELStripSizeY[n]);
    h18[1][n]->GetXaxis()->SetTitle("Y_{Cog}-Y_{Dig}");
    h18[1][n]->GetYaxis()->SetTitle("Number of events");
    sprintf(TELRoot,"Cog_Y_%u",n+1);
    h18[1][n]->SetTitle(TELRoot);
  }

  //Histograms for correlation
  TH2F* h19[2][glob->TELNumberOfPlanes - 1];
  TCanvas* canv19 = new TCanvas("c9","",600,600);
  canv19->SetFillColor(0);
  canv19->Divide(glob->TELNumberOfPlanes - 1,2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){

    sprintf(TELRoot,"WinCorrX_%u",n+1);
    h19[0][n] = new TH2F(TELRoot,"",glob->TELNumberOfStripsX[n],0,glob->TELNumberOfStripsX[n] * glob->TELStripSizeX[n],glob->TELNumberOfStripsX[n+1],0,glob->TELNumberOfStripsX[n+1] * glob->TELStripSizeX[n+1]); 
    sprintf(TELRoot,"X telescope %d",n+1);
    h19[0][n]->GetXaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"X telescope %d",n+2);
    h19[0][n]->GetYaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Correlation X_%u",n+1);
    h19[0][n]->SetTitle(TELRoot);

    sprintf(TELRoot,"WinCorrY_%u",n+1);
    h19[1][n] = new TH2F(TELRoot,"",glob->TELNumberOfStripsY[n],0,glob->TELNumberOfStripsY[n] * glob->TELStripSizeY[n],glob->TELNumberOfStripsY[n+1],0,glob->TELNumberOfStripsY[n+1] * glob->TELStripSizeY[n+1]); 
    sprintf(TELRoot,"Y telescope %d",n+1);
    h19[1][n]->GetXaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Y telescope %d",n+2);
    h19[1][n]->GetYaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Correlation Y_%u",n+1);
    h19[1][n]->SetTitle(TELRoot);

  }

  for(Int_t i = 0; i < ntuple->GetEntries(); i++){

    ntuple->GetEntry(i);

    for(Int_t i = 0; i < ADCXSize; i++){
      SeedX[i] = false;
      HitBufferX[i] = 0.0;
      CogBufferX[i] = 0.0;
    }
    for(Int_t i = 0; i < ADCYSize; i++){
      SeedY[i] = false;
      HitBufferY[i] = 0.0;
      CogBufferY[i] = 0.0;
    }

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

      for(Int_t x = ClustSizeX/2; x <= glob->TELNumberOfStripsX[n] - ClustSizeX/2; x++){

	if(PlaneXADC[x + ShiftX[n]] > glob->TELSeedCutX[n]){

	  SeedX[x + ShiftX[n]] = true;
	  
	  for(Int_t tmp_pointer = -ClustSizeX/2; tmp_pointer <= ClustSizeX/2; tmp_pointer++){
	    if(PlaneXADC[x + ShiftX[n]] < PlaneXADC[x + tmp_pointer + ShiftX[n]]){
	      SeedX[x + ShiftX[n]] = false;
	      continue;
	    }
	  }//tmp_pointer

	}//if

      }//x

      NumberOfHitsX = 0;

      for(Int_t x = 0; x < glob->TELNumberOfStripsX[n]; x++){

	if(SeedX[x + ShiftX[n]]){

	  HitBufferX[NumberOfHitsX + ShiftX[n]] = (Float_t)x * glob->TELStripSizeX[n];

	  ClustCharge = 0.0;
	  
	  for(Int_t tmp_pointer = -ClustSizeX/2; tmp_pointer <= ClustSizeX/2; tmp_pointer++){
	    ClustCharge += (Float_t)PlaneXADC[x + tmp_pointer + ShiftX[n]]; 
	    CogBufferX[NumberOfHitsX + ShiftX[n]] += (Float_t)(x + tmp_pointer)*PlaneXADC[x + tmp_pointer + ShiftX[n]];
	  }
	  
	  if(ClustCharge > 0.0){
	    h14[0][n]->Fill(ClustCharge);
	    CogBufferX[NumberOfHitsX + ShiftX[n]] = glob->TELStripSizeX[n] * CogBufferX[NumberOfHitsX + ShiftX[n]]/ClustCharge;
	  }

	  NumberOfHitsX++;

	}//if

      }//x

      h15[0][n]->Fill(NumberOfHitsX);

      for(Int_t y = ClustSizeY/2; y <= glob->TELNumberOfStripsY[n] - ClustSizeY/2; y++){

	if(PlaneYADC[y + ShiftY[n]] > glob->TELSeedCutY[n]){

	  SeedY[y + ShiftY[n]] = true;

	  for(Int_t tmp_pointer = -ClustSizeY/2; tmp_pointer <= ClustSizeY/2; tmp_pointer++){
	    if(PlaneYADC[y + ShiftY[n]] < PlaneYADC[y + tmp_pointer + ShiftY[n]]){
	      SeedY[y + ShiftY[n]] = false;
	      continue;
	    }
	  }//tmp_pointer

	}//if

      }//y

      NumberOfHitsY = 0;

      for(Int_t y = 0; y < glob->TELNumberOfStripsY[n]; y++){

	if(SeedY[y + ShiftY[n]]){
	  
	  HitBufferY[NumberOfHitsY + ShiftY[n]] = (Float_t)y * glob->TELStripSizeY[n];

	  ClustCharge = 0.0;

	  for(Int_t tmp_pointer = -ClustSizeY/2; tmp_pointer <= ClustSizeY/2; tmp_pointer++){
	    ClustCharge += (Float_t)PlaneYADC[y + tmp_pointer + ShiftY[n]]; 
	    CogBufferY[NumberOfHitsY + ShiftY[n]] += (Float_t)(y + tmp_pointer)*PlaneYADC[y + tmp_pointer + ShiftY[n]];
	  }
	  
	  if(ClustCharge > 0.0){
	    h14[1][n]->Fill(ClustCharge);
	    CogBufferY[NumberOfHitsY + ShiftY[n]] = glob->TELStripSizeY[n] * CogBufferY[NumberOfHitsY + ShiftY[n]]/ClustCharge;
	  }

	  NumberOfHitsY++;

	}//if

      }//y

      h15[1][n]->Fill(NumberOfHitsY);

      for(Int_t x = 0; x < NumberOfHitsX; x++){
	for(Int_t y = 0; y < NumberOfHitsY; y++){
	  h16[n]->Fill(HitBufferX[x + ShiftX[n]],HitBufferY[y + ShiftY[n]],1.0);
	  //h16[n]->Fill(CogBufferX[x + ShiftX[n]],CogBufferY[y + ShiftY[n]],1.0);
	}
      }

    }//n

    Int_t while1;
    Int_t while2;
    Float_t Distance;
    Float_t TmpDigX;
    Float_t TmpDigY;
    Float_t TmpCogX;
    Float_t TmpCogY;
    
    for(while1 = 0; while1 < glob->TELMaxHit; while1++){
      
      for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
	DigX[n][while1] = -1.0;
	CogX[n][while1] = -1.0;
	DigY[n][while1] = -1.0;
	CogY[n][while1] = -1.0;
      }

    }

    for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){

      Distance = glob->TELCorrIntervalX[n];

      while1 = 0;
      while2 = 0;
      
      while(CogBufferX[while1 + ShiftX[n]] > 0.0 && while1 < glob->TELMaxHit){ 


	if(n == 0){

	  DigX[n][while1] = HitBufferX[while1 + ShiftX[n]];
	  CogX[n][while1] = CogBufferX[while1 + ShiftX[n]];
	  h18[0][n]->Fill(CogX[n][while1] - DigX[n][while1]);

	}

	TmpDigX = -1.0;
	TmpCogX = -1.0;

	while(CogBufferX[while2 + ShiftX[n + 1]] > 0.0){
	  
	  h17[0][n]->Fill(DigX[n][while1],HitBufferX[while2 + ShiftX[n + 1]]);	  

	  if(TMath::Abs(glob->TELCorrealationSlopeX[n] * DigX[n][while1] + glob->TELCorrealationConstX[n] - HitBufferX[while2 + ShiftX[n + 1]]) < Distance){
	    Distance = TMath::Abs(glob->TELCorrealationSlopeX[n] * DigX[n][while1] + glob->TELCorrealationConstX[n] - HitBufferX[while2 + ShiftX[n + 1]]);
	    TmpDigX = HitBufferX[while2 + ShiftX[n + 1]];
	    TmpCogX = CogBufferX[while2 + ShiftX[n + 1]];
	  }

	  while2++;
	}//while2
	
	DigX[n + 1][while1] = TmpDigX;
	CogX[n + 1][while1] = TmpCogX;
	if(CogX[n + 1][while1] > 0.0) h18[0][n + 1]->Fill(CogX[n + 1][while1] - DigX[n + 1][while1]);
	
	while1++;
      }//while1

      if(while1 >= glob->TELMaxHit) cout << " WORNING: \n    Tom many hits in the x plane nr " << n  << " in event " << (Int_t)i << endl;

      Distance = glob->TELCorrIntervalY[n];

      while1 = 0;
      while2 = 0;
      
      while(CogBufferY[while1 + ShiftY[n]] > 0.0 && while1 < glob->TELMaxHit){ 

	if(n == 0){

	  DigY[n][while1] = HitBufferY[while1 + ShiftY[n]];
	  CogY[n][while1] = CogBufferY[while1 + ShiftY[n]];
	  h18[1][n]->Fill(CogY[n][while1] - DigY[n][while1]);

	}

	TmpDigY = -1.0;
	TmpCogY = -1.0;

	while(CogBufferY[while2 + ShiftY[n + 1]] > 0.0){
	  
	  h17[1][n]->Fill(DigY[n][while1],HitBufferY[while2 + ShiftY[n + 1]]);
		  
          if(TMath::Abs(glob->TELCorrealationSlopeY[n] * DigY[n][while1] + glob->TELCorrealationConstY[n] - HitBufferY[while2 + ShiftY[n + 1]]) < Distance){
	    Distance = TMath::Abs(glob->TELCorrealationSlopeY[n] * DigY[n][while1] + glob->TELCorrealationConstY[n] - HitBufferY[while2 + ShiftY[n + 1]]);
	    TmpDigY = HitBufferY[while2 + ShiftY[n + 1]];
	    TmpCogY = CogBufferY[while2 + ShiftY[n + 1]];
          }

	  while2++;
	}//while2

	DigY[n + 1][while1] = TmpDigY;
	CogY[n + 1][while1] = TmpCogY;
	if(CogY[n + 1][while1] > 0.0) h18[1][n + 1]->Fill(CogY[n + 1][while1] - DigY[n + 1][while1]);

	while1++;
      }//while1      

      if(while1 >= glob->TELMaxHit) cout << " WORNING: \n    Tom many hits in the y plane nr " << n  << " in event " << (Int_t)i << endl;

    }//n

    for(while1 = 0; while1 < glob->TELNumberOfPlanes - 1; while1++){

      VecFill = true;

      for(Int_t n = while1; n <= while1 + 1; n++){
	
	if(DigX[n][0] < 0.0  || DigX[n][1] > 0.0 || DigY[n][0] < 0.0  || DigY[n][1] > 0.0){
	  VecFill = false;
	  break;
	}
	
      }//n
      
      if(VecFill){
	
	for(Int_t n = while1; n <= while1 + 1; n++){
	  
	  VecDigX[while1][n].push_back((Double_t)DigX[n][0]);
	  VecDigY[while1][n].push_back((Double_t)DigY[n][0]);
	  VecCogX[while1][n].push_back((Double_t)CogX[n][0]);
	  VecCogY[while1][n].push_back((Double_t)CogY[n][0]);
	  
	}//n
	
      }//if

    }//while1

  }//i

  for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){

    for(Int_t i = 0; i < 6; i++) APara[i + n*6] = 0.0; 

    if(n > 0){

      for(Int_t i = 0; i < VecDigX[n][n].size(); i++) VecDigX[n][n][i] = TransX(VecDigX[n][n][i],VecDigY[n][n][i],&APara[(n-1)*6]); 
      for(Int_t i = 0; i < VecDigY[n][n].size(); i++) VecDigY[n][n][i] = TransY(VecDigX[n][n][i],VecDigY[n][n][i],&APara[(n-1)*6]); 
      for(Int_t i = 0; i < VecCogX[n][n].size(); i++) VecCogX[n][n][i] = TransX(VecCogX[n][n][i],VecCogY[n][n][i],&APara[(n-1)*6]); 
      for(Int_t i = 0; i < VecCogY[n][n].size(); i++) VecCogY[n][n][i] = TransY(VecCogX[n][n][i],VecCogY[n][n][i],&APara[(n-1)*6]); 

    }

    /*
    PreAlignment = new TELAlign(&VecDigX[n][n],&VecDigX[n][n+1],&VecDigY[n][n],&VecDigY[n][n+1]);
    PreAlignment->Align(&APara[n*6]);

    for(Int_t i = 0; i < VecDigX[n][n+1].size(); i++) VecDigX[n][n+1][i] = TransX(VecDigX[n][n+1][i],VecDigY[n][n+1][i],&APara[n*6]); 
    for(Int_t i = 0; i < VecDigY[n][n+1].size(); i++) VecDigY[n][n+1][i] = TransY(VecDigX[n][n+1][i],VecDigY[n][n+1][i],&APara[n*6]); 

    for(Int_t i = 0; i < VecDigX[n][n].size(); i++) h19[0][n]->Fill(VecDigX[n][n][i],VecDigX[n][n+1][i]);	  
    for(Int_t i = 0; i < VecDigY[n][n].size(); i++) h19[1][n]->Fill(VecDigY[n][n][i],VecDigY[n][n+1][i]);
    */

    PreAlignment = new TELAlign(&VecCogX[n][n],&VecCogX[n][n+1],&VecCogY[n][n],&VecCogY[n][n+1]);
    PreAlignment->Align(&APara[n*6]);

    for(Int_t i = 0; i < VecCogX[n][n+1].size(); i++) VecCogX[n][n+1][i] = TransX(VecCogX[n][n+1][i],VecCogY[n][n+1][i],&APara[n*6]); 
    for(Int_t i = 0; i < VecCogY[n][n+1].size(); i++) VecCogY[n][n+1][i] = TransY(VecCogX[n][n+1][i],VecCogY[n][n+1][i],&APara[n*6]); 

    for(Int_t i = 0; i < VecCogX[n][n].size(); i++) h19[0][n]->Fill(VecCogX[n][n][i],VecCogX[n][n+1][i]);	  
    for(Int_t i = 0; i < VecCogY[n][n].size(); i++) h19[1][n]->Fill(VecCogY[n][n][i],VecCogY[n][n+1][i]);

  }

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    canv14->cd(n+1);
    h14[0][n]->Draw();
    h14[0][n]->Fit("landau");
    canv14->cd(n + 1  + glob->TELNumberOfPlanes);
    h14[1][n]->Draw();
    h14[1][n]->Fit("landau");

    canv15->cd(n+1);
    h15[0][n]->Draw();
    canv15->cd(n + 1  + glob->TELNumberOfPlanes);
    h15[1][n]->Draw();

    canv16->cd(n+1);
    h16[n]->Draw("COLZ");

    canv18->cd(n+1);
    h18[0][n]->Draw();
    canv18->cd(n + 1  + glob->TELNumberOfPlanes);
    h18[1][n]->Draw();

  }

  /*
  TCanvas* PHDCanv1 = new TCanvas("PHDcanv1","",800,800);
  PHDCanv1->SetFillColor(0);

  PHDCanv1->cd();
  h14[0][0]->Draw();
  h14[0][0]->Fit("landau");

  TCanvas* PHDCanv2 = new TCanvas("PHDcanv2","",800,800);
  PHDCanv2->SetFillColor(0);

  PHDCanv2->cd();
  h14[1][0]->Draw();
  h14[1][0]->Fit("landau");
  */

  TF1* CorrelatonX[glob->TELNumberOfPlanes - 1];
  TF1* CorrelatonY[glob->TELNumberOfPlanes - 1];

  for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){

    sprintf(TELRoot,"f1X%u",n+1);
    CorrelatonX[n] = new TF1(TELRoot,"[0]*x + [1]",0,glob->TELNumberOfStripsX[n]*glob->TELStripSizeX[n]);
    CorrelatonX[n]->SetLineWidth(1);
    CorrelatonX[n]->SetLineColor(2); 
    CorrelatonX[n]->SetParameter(0,glob->TELCorrealationSlopeX[n]);
    CorrelatonX[n]->SetParameter(1,glob->TELCorrealationConstX[n]);
    canv17->cd(n+1);
    h17[0][n]->Draw("COLZ");
    CorrelatonX[n]->Draw("SAME");

    sprintf(TELRoot,"f1Y%u",n+1);
    CorrelatonY[n] = new TF1(TELRoot,"[0]*x + [1]",0,glob->TELNumberOfStripsY[n]*glob->TELStripSizeY[n]);
    CorrelatonY[n]->SetLineWidth(1);
    CorrelatonY[n]->SetLineColor(2);
    CorrelatonY[n]->SetParameter(0,glob->TELCorrealationSlopeY[n]);
    CorrelatonY[n]->SetParameter(1,glob->TELCorrealationConstY[n]);
    canv17->cd(n + glob->TELNumberOfPlanes);
    h17[1][n]->Draw("COLZ");
    CorrelatonY[n]->Draw("SAME");

    canv19->cd(n+1);
    canv19->GetPad(n+1)->SetGridx();
    canv19->GetPad(n+1)->SetGridy();
    h19[0][n]->Draw("COLZ");
    canv19->cd(n + glob->TELNumberOfPlanes);
    canv19->GetPad(n + glob->TELNumberOfPlanes)->SetGridx();
    canv19->GetPad(n + glob->TELNumberOfPlanes)->SetGridy();
    h19[1][n]->Draw("COLZ");

  }

  delete [] HitBufferX;
  delete [] HitBufferY;
  delete [] CogBufferX;
  delete [] CogBufferY;
  delete [] SeedX;
  delete [] SeedY;

}//FindClusters

void TELFile::HitNtupl(Int_t ClustSizeX, Int_t ClustSizeY){
  
  //Histograms for correlation
  TH2F* h20[2][glob->TELNumberOfPlanes - 1];
  TCanvas* canv20 = new TCanvas("c20","",600,600);
  canv20->SetFillColor(0);
  canv20->Divide(glob->TELNumberOfPlanes - 1,2);
  
  for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){
    
    sprintf(TELRoot,"WinCorrXHitNtupl_%u",n+1);
    h20[0][n] = new TH2F(TELRoot,"",glob->TELNumberOfStripsX[n],0,glob->TELNumberOfStripsX[n] * glob->TELStripSizeX[n],glob->TELNumberOfStripsX[n+1],0,glob->TELNumberOfStripsX[n+1] * glob->TELStripSizeX[n+1]); 
    sprintf(TELRoot,"X telescope %d",n+1);
    h20[0][n]->GetXaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"X telescope %d",n+2);
    h20[0][n]->GetYaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Correlation X_%u",n+1);
    h20[0][n]->SetTitle(TELRoot);

    sprintf(TELRoot,"WinCorrYHitNtupl_%u",n+1);
    h20[1][n] = new TH2F(TELRoot,"",glob->TELNumberOfStripsY[n],0,glob->TELNumberOfStripsY[n] * glob->TELStripSizeY[n],glob->TELNumberOfStripsY[n+1],0,glob->TELNumberOfStripsY[n+1] * glob->TELStripSizeY[n+1]); 
    sprintf(TELRoot,"Y telescope %d",n+1);
    h20[1][n]->GetXaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Y telescope %d",n+2);
    h20[1][n]->GetYaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Correlation Y_%u",n+1);
    h20[1][n]->SetTitle(TELRoot);

  }

  Float_t DigX[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Float_t DigY[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Float_t CogX[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Float_t CogY[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Double_t TrigTime[2];

  sprintf(TELRoot,"ntuple/tele/TELHits%u.root",RunNumber);

  TFile *trackfile;
  trackfile = new TFile(TELRoot,"RECREATE");
  trackfile->cd();

  TTree *Tracks = new TTree("TELHits","filename");

  Tracks->Branch("Triger",&TrigTime[0],"TrigTime[2]/D");
  sprintf(TELRoot,"DigX[%u][%u]",glob->TELNumberOfPlanes,glob->TELMaxHit);
  Tracks->Branch("DigX",&DigX[0][0],TELRoot);
  sprintf(TELRoot,"DigY[%u][%u]",glob->TELNumberOfPlanes,glob->TELMaxHit);
  Tracks->Branch("DigY",&DigY[0][0],TELRoot);
  sprintf(TELRoot,"CogX[%u][%u]",glob->TELNumberOfPlanes,glob->TELMaxHit);
  Tracks->Branch("CogX",&CogX[0][0],TELRoot);
  sprintf(TELRoot,"CogY[%u][%u]",glob->TELNumberOfPlanes,glob->TELMaxHit);
  Tracks->Branch("CogY",&CogY[0][0],TELRoot);

  vector<Double_t> VecDigX[glob->TELNumberOfPlanes];
  vector<Double_t> VecDigY[glob->TELNumberOfPlanes];
  vector<Double_t> VecCogX[glob->TELNumberOfPlanes];
  vector<Double_t> VecCogY[glob->TELNumberOfPlanes];
  Bool_t VecFill;
  
  Bool_t* SeedX = new Bool_t[ADCXSize];
  Bool_t* SeedY = new Bool_t[ADCYSize];
  
  Float_t* HitBufferX = new Float_t[ADCXSize];
  Float_t* HitBufferY = new Float_t[ADCYSize];
  Float_t* CogBufferX = new Float_t[ADCXSize];
  Float_t* CogBufferY = new Float_t[ADCYSize];
  
  sprintf(TELRoot,"ntuple/tele/TELData%u.root",RunNumber);
  
  TELData = new TFile(TELRoot);
  TTree* ntuple = (TTree*)TELData->Get("TELData");
  
  ntuple->SetBranchAddress("Trigger",&Trigger[0]);
  ntuple->SetBranchAddress("PlaneX",&PlaneXADC[0]);
  ntuple->SetBranchAddress("PlaneY",&PlaneYADC[0]);  

  Float_t ClustCharge;  

  Int_t NumberOfHitsX;
  Int_t NumberOfHitsY;

  for(Int_t i = 0; i < ntuple->GetEntries(); i++){
    
    ntuple->GetEntry(i);
    
    for(Int_t i = 0; i < ADCXSize; i++){
      SeedX[i] = false;
      HitBufferX[i] = 0.0;
      CogBufferX[i] = 0.0;
    }
    for(Int_t i = 0; i < ADCYSize; i++){
      SeedY[i] = false;
      HitBufferY[i] = 0.0;
      CogBufferY[i] = 0.0;
    }

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

      for(Int_t x = ClustSizeX/2; x <= glob->TELNumberOfStripsX[n] - ClustSizeX/2; x++){

	if(PlaneXADC[x + ShiftX[n]] > glob->TELSeedCutX[n]){

	  SeedX[x + ShiftX[n]] = true;

	  for(Int_t tmp_pointer = -ClustSizeX/2; tmp_pointer <= ClustSizeX/2; tmp_pointer++){
	    if(PlaneXADC[x + ShiftX[n]] < PlaneXADC[x + tmp_pointer + ShiftX[n]]){
	      SeedX[x + ShiftX[n]] = false;
	      continue;
	    }
	  }//tmp_pointer

	}//if

      }//x

      NumberOfHitsX = 0;

      for(Int_t x = 0; x < glob->TELNumberOfStripsX[n]; x++){

	if(SeedX[x + ShiftX[n]]){

	  HitBufferX[NumberOfHitsX + ShiftX[n]] = (Float_t)x * glob->TELStripSizeX[n];

	  ClustCharge = 0.0;
	  
	  for(Int_t tmp_pointer = -ClustSizeX/2; tmp_pointer <= ClustSizeX/2; tmp_pointer++){
	    ClustCharge += (Float_t)PlaneXADC[x + tmp_pointer + ShiftX[n]]; 
	    CogBufferX[NumberOfHitsX + ShiftX[n]] += (Float_t)(x + tmp_pointer)*PlaneXADC[x + tmp_pointer + ShiftX[n]];
	  }
	  
	  if(ClustCharge > 0.0){
	    CogBufferX[NumberOfHitsX + ShiftX[n]] = glob->TELStripSizeX[n] * CogBufferX[NumberOfHitsX + ShiftX[n]]/ClustCharge;
	  }

	  NumberOfHitsX++;

	}//if

      }//x

      for(Int_t y = ClustSizeY/2; y <= glob->TELNumberOfStripsY[n] - ClustSizeY/2; y++){

	if(PlaneYADC[y + ShiftY[n]] > glob->TELSeedCutY[n]){

	  SeedY[y + ShiftY[n]] = true;

	  for(Int_t tmp_pointer = -ClustSizeY/2; tmp_pointer <= ClustSizeY/2; tmp_pointer++){
	    if(PlaneYADC[y + ShiftY[n]] < PlaneYADC[y + tmp_pointer + ShiftY[n]]){
	      SeedY[y + ShiftY[n]] = false;
	      continue;
	    }
	  }//tmp_pointer

	}//if

      }//y

      NumberOfHitsY = 0;

      for(Int_t y = 0; y < glob->TELNumberOfStripsY[n]; y++){

	if(SeedY[y + ShiftY[n]]){
	  
	  HitBufferY[NumberOfHitsY + ShiftY[n]] = (Float_t)y * glob->TELStripSizeY[n];

	  ClustCharge = 0.0;

	  for(Int_t tmp_pointer = -ClustSizeY/2; tmp_pointer <= ClustSizeY/2; tmp_pointer++){
	    ClustCharge += (Float_t)PlaneYADC[y + tmp_pointer + ShiftY[n]]; 
	    CogBufferY[NumberOfHitsY + ShiftY[n]] += (Float_t)(y + tmp_pointer)*PlaneYADC[y + tmp_pointer + ShiftY[n]];
	  }
	  
	  if(ClustCharge > 0.0){
	    CogBufferY[NumberOfHitsY + ShiftY[n]] = glob->TELStripSizeY[n] * CogBufferY[NumberOfHitsY + ShiftY[n]]/ClustCharge;
	  }

	  NumberOfHitsY++;

	}//if

      }//y

    }//n

    Int_t while1;
    Int_t while2;
    Float_t Distance;
    
    for(while1 = 0; while1 < glob->TELMaxHit; while1++){
      
      for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
	DigX[n][while1] = -1.0;
	CogX[n][while1] = -1.0;
	DigY[n][while1] = -1.0;
	CogY[n][while1] = -1.0;
      }

    }

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

      while1 = 0;
      while2 = 0;
	
      while(CogBufferX[while1 + ShiftX[n]] > 0.0 && while1 < glob->TELMaxHit){ 

	if(n == 0){
	  DigX[n][while1] = ((Double_t)HitBufferX[while1 + ShiftX[n]]);
	  CogX[n][while1] = ((Double_t)CogBufferX[while1 + ShiftX[n]]);
	}

	while(CogBufferY[while2 + ShiftY[n]] > 0.0 && while2 < glob->TELMaxHit){
	  
	  if(n == 0){
	    DigY[n][while2] = ((Double_t)HitBufferY[while2 + ShiftY[n]]);
	    CogY[n][while2] = ((Double_t)CogBufferY[while2 + ShiftY[n]]);
	  }
	  if(n > 0){
	    VecDigX[n].push_back(TransX((Double_t)HitBufferX[while1 + ShiftX[n]],(Double_t)HitBufferY[while2 + ShiftY[n]],&APara[(n-1)*6]));
	    VecDigY[n].push_back(TransY((Double_t)HitBufferX[while1 + ShiftX[n]],(Double_t)HitBufferY[while2 + ShiftY[n]],&APara[(n-1)*6]));
	    VecCogX[n].push_back(TransX((Double_t)CogBufferX[while1 + ShiftX[n]],(Double_t)CogBufferY[while2 + ShiftY[n]],&APara[(n-1)*6]));
	    VecCogY[n].push_back(TransY((Double_t)CogBufferX[while1 + ShiftX[n]],(Double_t)CogBufferY[while2 + ShiftY[n]],&APara[(n-1)*6])); 
	  }

	  while2++;

	}//while2

	while1++;

      }//while1

    }//n

    for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){

      while1 = 0;
      while(DigX[n][while1] > 0.0 && while1 < glob->TELMaxHit){

	Distance = glob->TELCorrIntervalX[n];

	for(while2 = 0; while2 < VecDigX[n+1].size(); while2++){    

	  if(TMath::Abs(DigX[n][while1] - VecDigX[n+1][while2]) < Distance){
	    Distance = TMath::Abs(DigX[n][while1] - VecDigX[n+1][while2]);
	    DigX[n+1][while1] = VecDigX[n+1][while2];
	    CogX[n+1][while1] = VecCogX[n+1][while2];
	  }

	}//while2

	while1++;

      }//while1

      while1 = 0;
      while(DigY[n][while1] > 0.0 && while1 < glob->TELMaxHit){

	Distance = glob->TELCorrIntervalY[n];

	for(while2 = 0; while2 < VecDigY[n+1].size(); while2++){

	  if(TMath::Abs(DigY[n][while1] - VecDigY[n+1][while2]) < Distance){
	    Distance = TMath::Abs(DigY[n][while1] - VecDigY[n+1][while2]);
	    DigY[n+1][while1] = VecDigY[n+1][while2];
	    CogY[n+1][while1] = VecCogY[n+1][while2];
	  }

	}//while2
	
	while1++;
	
      }//while1
      
    }//n

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
      
      VecDigX[n].clear();
      VecDigY[n].clear();
      VecCogX[n].clear();
      VecCogY[n].clear();

    }

    for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){  

      while1 = 0;
      while(DigX[n][while1] > 0.0){
	h20[0][n]->Fill(DigX[n][while1],DigX[n+1][while1]);	
	//h20[0][n]->Fill(CogX[n][while1],CogX[n+1][while1]);	
	while1++;
      }

      while1 = 0;
      while(DigY[n][while1] > 0.0){
	h20[1][n]->Fill(DigY[n][while1],DigY[n+1][while1]);	
	//h20[1][n]->Fill(CogY[n][while1],CogY[n+1][while1]);	
	while1++;
      }

    }//n

    TrigTime[0] = Trigger[0];
    TrigTime[1] = Trigger[1];

    Tracks->Fill();

  }//i

  trackfile->cd();
  Tracks->Write();
  trackfile->Close();

  for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){

    canv20->cd(n+1);
    canv20->GetPad(n+1)->SetGridx();
    canv20->GetPad(n+1)->SetGridy();
    h20[0][n]->Draw("COLZ");
    canv20->cd(n + glob->TELNumberOfPlanes);
    canv20->GetPad(n + glob->TELNumberOfPlanes)->SetGridx();
    canv20->GetPad(n + glob->TELNumberOfPlanes)->SetGridy();
    h20[1][n]->Draw("COLZ");

  }
  
}//HitNtupl

void TELFile::HitNtuplWithEta(Int_t ClustSizeX, Int_t ClustSizeY, Int_t Binning){
  
  //Histograms for signals
  TH1F* h10[2][glob->TELNumberOfPlanes];
  TCanvas* canv10 = new TCanvas("c10","",1200,800);
  canv10->SetFillColor(0);
  canv10->Divide(glob->TELNumberOfPlanes,2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
    sprintf(TELRoot,"Eta_X_%u",n+1);
    h10[0][n] = new TH1F(TELRoot,"",Binning,0,1); 
    h10[0][n]->GetXaxis()->SetTitle("#eta");
    h10[0][n]->GetYaxis()->SetTitle("Number of events");
    h10[0][n]->GetYaxis()->SetTitle("");
    sprintf(TELRoot,"Telescope %u (plane X)",n+1);
    h10[0][n]->SetTitle(TELRoot);
    h10[0][n]->SetLineWidth(3);
    h10[0][n]->SetLineColor(4);

    sprintf(TELRoot,"Eta_Y_%u",n+1);
    h10[1][n] = new TH1F(TELRoot,"",Binning,0,1);
    h10[1][n]->GetXaxis()->SetTitle("#eta");
    h10[1][n]->GetYaxis()->SetTitle("Number of events");
    h10[1][n]->GetYaxis()->SetTitle("");
    sprintf(TELRoot,"Telescope %u (plane Y)",n+1);
    h10[1][n]->SetTitle(TELRoot);
    h10[1][n]->SetLineWidth(3);
    h10[1][n]->SetLineColor(4);
  }//n

  //Histograms for signals
  TH1F* h11[2][glob->TELNumberOfPlanes];
  TCanvas* canv11 = new TCanvas("c11","",1200,800);
  canv11->SetFillColor(0);
  canv11->Divide(glob->TELNumberOfPlanes,2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    sprintf(TELRoot,"EtaDist_X_%u",n+1);
    h11[0][n] = new TH1F(TELRoot,"",Binning,0,1); 
    h11[0][n]->GetXaxis()->SetTitle("#eta");
    h11[0][n]->GetYaxis()->SetTitle("f(#eta)");
    sprintf(TELRoot,"Telescope %u (plane X)",n+1);
    h11[0][n]->SetTitle(TELRoot);
    h11[0][n]->SetLineWidth(3);
    h11[0][n]->SetLineColor(4);
    h11[0][n]->SetStats(0);

    sprintf(TELRoot,"EtaDist_Y_%u",n+1);
    h11[1][n] = new TH1F(TELRoot,"",Binning,0,1);
    h11[1][n]->GetXaxis()->SetTitle("#eta");
    h11[1][n]->GetYaxis()->SetTitle("f(#eta)");
    sprintf(TELRoot,"Telescope %u (plane Y)",n+1);
    h11[1][n]->SetTitle(TELRoot);
    h11[1][n]->SetLineWidth(3);
    h11[1][n]->SetLineColor(4);
    h11[1][n]->SetStats(0);

  }//n

  Float_t EtaCorrectionX[glob->TELNumberOfPlanes][Binning];
  Float_t EtaCorrectionY[glob->TELNumberOfPlanes][Binning];
  
  //Histograms for correlation
  TH2F* h20[2][glob->TELNumberOfPlanes - 1];
  TCanvas* canv20 = new TCanvas("c20","",600,600);
  canv20->SetFillColor(0);
  canv20->Divide(glob->TELNumberOfPlanes - 1,2);
  
  for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){
    
    sprintf(TELRoot,"WinCorrXHitNtupl_%u",n+1);
    h20[0][n] = new TH2F(TELRoot,"",glob->TELNumberOfStripsX[n],0,glob->TELNumberOfStripsX[n] * glob->TELStripSizeX[n],glob->TELNumberOfStripsX[n+1],0,glob->TELNumberOfStripsX[n+1] * glob->TELStripSizeX[n+1]); 
    sprintf(TELRoot,"X telescope %d",n+1);
    h20[0][n]->GetXaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"X telescope %d",n+2);
    h20[0][n]->GetYaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Correlation X_%u",n+1);
    h20[0][n]->SetTitle(TELRoot);
    
    sprintf(TELRoot,"WinCorrYHitNtupl_%u",n+1);
    h20[1][n] = new TH2F(TELRoot,"",glob->TELNumberOfStripsY[n],0,glob->TELNumberOfStripsY[n] * glob->TELStripSizeY[n],glob->TELNumberOfStripsY[n+1],0,glob->TELNumberOfStripsY[n+1] * glob->TELStripSizeY[n+1]); 
    sprintf(TELRoot,"Y telescope %d",n+1);
    h20[1][n]->GetXaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Y telescope %d",n+2);
    h20[1][n]->GetYaxis()->SetTitle(TELRoot);
    sprintf(TELRoot,"Correlation Y_%u",n+1);
    h20[1][n]->SetTitle(TELRoot);
    
  }//n
  
  Float_t DigX[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Float_t DigY[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Float_t CogX[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Float_t CogY[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Double_t TrigTime[2];

  Float_t DigEtaX[glob->TELNumberOfPlanes][glob->TELMaxHit];
  Float_t DigEtaY[glob->TELNumberOfPlanes][glob->TELMaxHit];

  Double_t S_L, S_R, X_L, Y_L;
  Double_t Eta;

  sprintf(TELRoot,"ntuple/tele/TELHits%u.root",RunNumber);

  TFile *trackfile;
  trackfile = new TFile(TELRoot,"RECREATE");
  trackfile->cd();

  TTree *Tracks = new TTree("TELHits","filename");

  Tracks->Branch("Triger",&TrigTime[0],"TrigTime[2]/D");
  sprintf(TELRoot,"DigX[%u][%u]",glob->TELNumberOfPlanes,glob->TELMaxHit);
  Tracks->Branch("DigX",&DigX[0][0],TELRoot);
  sprintf(TELRoot,"DigY[%u][%u]",glob->TELNumberOfPlanes,glob->TELMaxHit);
  Tracks->Branch("DigY",&DigY[0][0],TELRoot);
  sprintf(TELRoot,"CogX[%u][%u]",glob->TELNumberOfPlanes,glob->TELMaxHit);
  Tracks->Branch("CogX",&CogX[0][0],TELRoot);
  sprintf(TELRoot,"CogY[%u][%u]",glob->TELNumberOfPlanes,glob->TELMaxHit);
  Tracks->Branch("CogY",&CogY[0][0],TELRoot);

  vector<Double_t> VecDigX[glob->TELNumberOfPlanes];
  vector<Double_t> VecDigY[glob->TELNumberOfPlanes];
  vector<Double_t> VecCogX[glob->TELNumberOfPlanes];
  vector<Double_t> VecCogY[glob->TELNumberOfPlanes];

  Bool_t VecFill;
  
  Bool_t* SeedX = new Bool_t[ADCXSize];
  Bool_t* SeedY = new Bool_t[ADCYSize];
  
  Float_t* HitBufferX = new Float_t[ADCXSize];
  Float_t* HitBufferY = new Float_t[ADCYSize];
  Float_t* CogBufferX = new Float_t[ADCXSize];
  Float_t* CogBufferY = new Float_t[ADCYSize];
  
  sprintf(TELRoot,"ntuple/tele/TELData%u.root",RunNumber);
  
  TELData = new TFile(TELRoot);
  TTree* ntuple = (TTree*)TELData->Get("TELData");
  
  ntuple->SetBranchAddress("Trigger",&Trigger[0]);
  ntuple->SetBranchAddress("PlaneX",&PlaneXADC[0]);
  ntuple->SetBranchAddress("PlaneY",&PlaneYADC[0]);  

  Float_t ClustCharge;  

  Int_t NumberOfHitsX;
  Int_t NumberOfHitsY;

  for(Int_t i = 0; i < ntuple->GetEntries(); i++){
    
    ntuple->GetEntry(i);
    
    for(Int_t i = 0; i < ADCXSize; i++){
      SeedX[i] = false;
      HitBufferX[i] = 0.0;
      CogBufferX[i] = 0.0;
    }
    for(Int_t i = 0; i < ADCYSize; i++){
      SeedY[i] = false;
      HitBufferY[i] = 0.0;
      CogBufferY[i] = 0.0;
    }

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
      
      for(Int_t x = ClustSizeX/2; x <= glob->TELNumberOfStripsX[n] - ClustSizeX/2; x++){

	if(PlaneXADC[x + ShiftX[n]] > glob->TELSeedCutX[n]){
	  
	  SeedX[x + ShiftX[n]] = true;

	  for(Int_t tmp_pointer = -ClustSizeX/2; tmp_pointer <= ClustSizeX/2; tmp_pointer++){
	    if(PlaneXADC[x + ShiftX[n]] < PlaneXADC[x + tmp_pointer + ShiftX[n]]){
	      SeedX[x + ShiftX[n]] = false;
	      continue;
	    }
	  }//tmp_pointer
	  
	}//if
      
      }//x
      
      NumberOfHitsX = 0;

      for(Int_t x = 0; x < glob->TELNumberOfStripsX[n]; x++){
	
	if(SeedX[x + ShiftX[n]]){
	  
	  HitBufferX[NumberOfHitsX + ShiftX[n]] = (Float_t)x * glob->TELStripSizeX[n];
	  
	  NumberOfHitsX++;
	  
	}//if

      }//x

      for(Int_t y = ClustSizeY/2; y <= glob->TELNumberOfStripsY[n] - ClustSizeY/2; y++){

	if(PlaneYADC[y + ShiftY[n]] > glob->TELSeedCutY[n]){

	  SeedY[y + ShiftY[n]] = true;

	  for(Int_t tmp_pointer = -ClustSizeY/2; tmp_pointer <= ClustSizeY/2; tmp_pointer++){
	    if(PlaneYADC[y + ShiftY[n]] < PlaneYADC[y + tmp_pointer + ShiftY[n]]){
	      SeedY[y + ShiftY[n]] = false;
	      continue;
	    }
	  }//tmp_pointer

	}//if

      }//y

      NumberOfHitsY = 0;

      for(Int_t y = 0; y < glob->TELNumberOfStripsY[n]; y++){

	if(SeedY[y + ShiftY[n]]){
	  
	  HitBufferY[NumberOfHitsY + ShiftY[n]] = (Float_t)y * glob->TELStripSizeY[n];
	  
	  NumberOfHitsY++;
	  
	}//if
	
      }//y
      
    }//n
    
    Int_t while1;
    Int_t while2;
    Float_t Distance;
    
    for(while1 = 0; while1 < glob->TELMaxHit; while1++){
      
      for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
	DigX[n][while1] = -1.0;
	CogX[n][while1] = -1.0;
	DigY[n][while1] = -1.0;
	CogY[n][while1] = -1.0;
	DigEtaX[n][while1] = -1.0;
	DigEtaY[n][while1] = -1.0;
      }

    }

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

      while1 = 0;
      while2 = 0;
      
      while(HitBufferX[while1 + ShiftX[n]] > 0.0 && while1 < glob->TELMaxHit){
	DigEtaX[n][while1] = HitBufferX[while1 + ShiftX[n]];
	while1++;
      }
      while(HitBufferY[while2 + ShiftY[n]] > 0.0 && while2 < glob->TELMaxHit){
	DigEtaY[n][while2] = HitBufferY[while2 + ShiftY[n]];
	while2++;
      }

      while1 = 0;
      while(DigEtaX[n][while1] > 0.0){
	
	S_R = 0.0;
	S_L =0.0;

	if(PlaneXADC[(Int_t)(DigEtaX[n][while1]/glob->TELStripSizeX[n]) - 1 + ShiftX[n]] < PlaneXADC[(Int_t)(DigEtaX[n][while1]/glob->TELStripSizeX[n]) + 1 + ShiftX[n]]){
	  S_L = PlaneXADC[(Int_t)(DigEtaX[n][while1]/glob->TELStripSizeX[n]) + ShiftX[n]];
	  S_R = PlaneXADC[(Int_t)(DigEtaX[n][while1]/glob->TELStripSizeX[n]) + 1 + ShiftX[n]];
	}
	if(PlaneXADC[(Int_t)(DigEtaX[n][while1]/glob->TELStripSizeX[n]) - 1 + ShiftX[n]] > PlaneXADC[(Int_t)(DigEtaX[n][while1]/glob->TELStripSizeX[n]) + 1 + ShiftX[n]]){
	  S_L = PlaneXADC[(Int_t)(DigEtaX[n][while1]/glob->TELStripSizeX[n]) - 1 + ShiftX[n]];
	  S_R = PlaneXADC[(Int_t)(DigEtaX[n][while1]/glob->TELStripSizeX[n]) + ShiftX[n]];
	}
	
	if(S_L > 0.0 && S_R > 0.0){
	  Eta = S_R/(S_R + S_L);
	  h10[0][n]->Fill(Eta);
	}

	while1++;
	
      }//while1

      while1 = 0;
      while(DigEtaY[n][while1] > 0.0){
	
	S_R = 0.0;
	S_L =0.0;

	if(PlaneYADC[(Int_t)(DigEtaY[n][while1]/glob->TELStripSizeY[n]) - 1 + ShiftY[n]] < PlaneYADC[(Int_t)(DigEtaY[n][while1]/glob->TELStripSizeY[n]) + 1 + ShiftY[n]]){
	  S_L = PlaneYADC[(Int_t)(DigEtaY[n][while1]/glob->TELStripSizeY[n]) + ShiftY[n]];
	  S_R = PlaneYADC[(Int_t)(DigEtaY[n][while1]/glob->TELStripSizeY[n]) + 1 + ShiftY[n]];
	}
	if(PlaneYADC[(Int_t)(DigEtaY[n][while1]/glob->TELStripSizeY[n]) - 1 + ShiftY[n]] > PlaneYADC[(Int_t)(DigEtaY[n][while1]/glob->TELStripSizeY[n]) + 1 + ShiftY[n]]){
	  S_L = PlaneYADC[(Int_t)(DigEtaY[n][while1]/glob->TELStripSizeY[n]) - 1 + ShiftY[n]];
	  S_R = PlaneYADC[(Int_t)(DigEtaY[n][while1]/glob->TELStripSizeY[n]) + ShiftY[n]];
	}
	
	if(S_L > 0.0 && S_R > 0.0){
	  Eta = S_R/(S_R + S_L);
	  h10[1][n]->Fill(Eta);
	}

	while1++;
	
      }//while1

    }//n

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
      
      EtaCorrectionX[n][0] = (Float_t)h10[0][n]->GetBinContent(1);
      EtaCorrectionY[n][0] = (Float_t)h10[1][n]->GetBinContent(1);
      
      for(Int_t b = 1; b < Binning; b++){
	
	EtaCorrectionX[n][b] = (Float_t)h10[0][n]->GetBinContent(b+1) + EtaCorrectionX[n][b-1];
	EtaCorrectionY[n][b] = (Float_t)h10[1][n]->GetBinContent(b+1) + EtaCorrectionY[n][b-1];

      }//b
      
    }//n

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
      for(Int_t b = 0; b < Binning; b++){

	EtaCorrectionX[n][b] = EtaCorrectionX[n][b]/EtaCorrectionX[n][Binning - 1];
	EtaCorrectionY[n][b] = EtaCorrectionY[n][b]/EtaCorrectionY[n][Binning - 1];

	h11[0][n]->SetBinContent(b+1,EtaCorrectionX[n][b]);
	h11[1][n]->SetBinContent(b+1,EtaCorrectionY[n][b]);

      }//b
    }//n
    
  }//i

  for(Int_t i = 0; i < ntuple->GetEntries(); i++){
    
    ntuple->GetEntry(i);
    
    for(Int_t i = 0; i < ADCXSize; i++){
      SeedX[i] = false;
      HitBufferX[i] = 0.0;
      CogBufferX[i] = 0.0;
    }
    for(Int_t i = 0; i < ADCYSize; i++){
      SeedY[i] = false;
      HitBufferY[i] = 0.0;
      CogBufferY[i] = 0.0;
    }

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
      
      for(Int_t x = ClustSizeX/2; x <= glob->TELNumberOfStripsX[n] - ClustSizeX/2; x++){
	
	if(PlaneXADC[x + ShiftX[n]] > glob->TELSeedCutX[n]){
	  
	  SeedX[x + ShiftX[n]] = true;
	  
	  for(Int_t tmp_pointer = -ClustSizeX/2; tmp_pointer <= ClustSizeX/2; tmp_pointer++){
	    if(PlaneXADC[x + ShiftX[n]] < PlaneXADC[x + tmp_pointer + ShiftX[n]]){
	      SeedX[x + ShiftX[n]] = false;
	      continue;
	    }
	  }//tmp_pointer
	  
	}//if
	
      }//x
      
      NumberOfHitsX = 0;
      
      for(Int_t x = 0; x < glob->TELNumberOfStripsX[n]; x++){
	
	if(SeedX[x + ShiftX[n]]){
	  
	  HitBufferX[NumberOfHitsX + ShiftX[n]] = (Float_t)x * glob->TELStripSizeX[n];
	  
	  S_R = 0.0;
	  S_L =0.0;
	  
	  if(PlaneXADC[x - 1 + ShiftX[n]] < PlaneXADC[x + 1 + ShiftX[n]]){
	    S_L = PlaneXADC[x + ShiftX[n]];
	    S_R = PlaneXADC[x + 1 + ShiftX[n]];
	    X_L = x * glob->TELStripSizeX[n];
	  }
	  if(PlaneXADC[x - 1 + ShiftX[n]] > PlaneXADC[x + 1 + ShiftX[n]]){
	    S_L = PlaneXADC[x - 1 + ShiftX[n]];
	    S_R = PlaneXADC[x + ShiftX[n]];
	    X_L = (x - 1) * glob->TELStripSizeX[n];
	  }
	  
	  if(S_L > 0.0 && S_R > 0.0) Eta = S_R/(S_R + S_L);

	  CogBufferX[NumberOfHitsX + ShiftX[n]] = X_L +  glob->TELStripSizeX[n] * EtaCorrectionX[n][(Int_t)(Eta*(Binning - 1))];

	  NumberOfHitsX++;
	  
	}//if
	
      }//x

      for(Int_t y = ClustSizeY/2; y <= glob->TELNumberOfStripsY[n] - ClustSizeY/2; y++){
	
	if(PlaneYADC[y + ShiftY[n]] > glob->TELSeedCutY[n]){
	  
	  SeedY[y + ShiftY[n]] = true;

	  for(Int_t tmp_pointer = -ClustSizeY/2; tmp_pointer <= ClustSizeY/2; tmp_pointer++){
	    if(PlaneYADC[y + ShiftY[n]] < PlaneYADC[y + tmp_pointer + ShiftY[n]]){
	      SeedY[y + ShiftY[n]] = false;
	      continue;
	    }
	  }//tmp_pointer

	}//if

      }//y

      NumberOfHitsY = 0;

      for(Int_t y = 0; y < glob->TELNumberOfStripsY[n]; y++){
	
	if(SeedY[y + ShiftY[n]]){
	  
	  HitBufferY[NumberOfHitsY + ShiftY[n]] = (Float_t)y * glob->TELStripSizeY[n];
	  
	  S_R = 0.0;
	  S_L = 0.0;
	  
	  if(PlaneYADC[y - 1 + ShiftY[n]] < PlaneYADC[y + 1 + ShiftY[n]]){
	    S_L = PlaneYADC[y + ShiftY[n]];
	    S_R = PlaneYADC[y + 1 + ShiftY[n]];
	    Y_L = y * glob->TELStripSizeY[n];
	  }
	  if(PlaneYADC[y - 1 + ShiftY[n]] > PlaneYADC[y + 1 + ShiftY[n]]){
	    S_L = PlaneYADC[y - 1 + ShiftY[n]];
	    S_R = PlaneYADC[y + ShiftY[n]];
	    Y_L = (y - 1) * glob->TELStripSizeY[n];
	  }
	  
	  if(S_L > 0.0 && S_R > 0.0) Eta = S_R/(S_R + S_L);

	  CogBufferY[NumberOfHitsY + ShiftY[n]] = Y_L +  glob->TELStripSizeY[n] * EtaCorrectionY[n][(Int_t)(Eta*(Binning - 1))];
	  
	  NumberOfHitsY++;

	}//if

      }//y

    }//n

    Int_t while1;
    Int_t while2;
    Float_t Distance;
    
    for(while1 = 0; while1 < glob->TELMaxHit; while1++){
      
      for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
	DigX[n][while1] = -1.0;
	CogX[n][while1] = -1.0;
	DigY[n][while1] = -1.0;
	CogY[n][while1] = -1.0;
      }

    }

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

      while1 = 0;
      while2 = 0;
	
      while(CogBufferX[while1 + ShiftX[n]] > 0.0 && while1 < glob->TELMaxHit){ 

	if(n == 0){
	  DigX[n][while1] = ((Double_t)HitBufferX[while1 + ShiftX[n]]);
	  CogX[n][while1] = ((Double_t)CogBufferX[while1 + ShiftX[n]]);
	}

	while(CogBufferY[while2 + ShiftY[n]] > 0.0 && while2 < glob->TELMaxHit){
	  
	  if(n == 0){
	    DigY[n][while2] = ((Double_t)HitBufferY[while2 + ShiftY[n]]);
	    CogY[n][while2] = ((Double_t)CogBufferY[while2 + ShiftY[n]]);
	  }
	  if(n > 0){
	    VecDigX[n].push_back(TransX((Double_t)HitBufferX[while1 + ShiftX[n]],(Double_t)HitBufferY[while2 + ShiftY[n]],&APara[(n-1)*6]));
	    VecDigY[n].push_back(TransY((Double_t)HitBufferX[while1 + ShiftX[n]],(Double_t)HitBufferY[while2 + ShiftY[n]],&APara[(n-1)*6]));
	    VecCogX[n].push_back(TransX((Double_t)CogBufferX[while1 + ShiftX[n]],(Double_t)CogBufferY[while2 + ShiftY[n]],&APara[(n-1)*6]));
	    VecCogY[n].push_back(TransY((Double_t)CogBufferX[while1 + ShiftX[n]],(Double_t)CogBufferY[while2 + ShiftY[n]],&APara[(n-1)*6])); 
	  }

	  while2++;

	}//while2

	while1++;

      }//while1

    }//n

    for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){

      while1 = 0;
      while(DigX[n][while1] > 0.0 && while1 < glob->TELMaxHit){

	Distance = glob->TELCorrIntervalX[n];

	for(while2 = 0; while2 < VecDigX[n+1].size(); while2++){    

	  if(TMath::Abs(DigX[n][while1] - VecDigX[n+1][while2]) < Distance){
	    Distance = TMath::Abs(DigX[n][while1] - VecDigX[n+1][while2]);
	    DigX[n+1][while1] = VecDigX[n+1][while2];
	    CogX[n+1][while1] = VecCogX[n+1][while2];
	  }

	}//while2

	while1++;

      }//while1

      while1 = 0;
      while(DigY[n][while1] > 0.0 && while1 < glob->TELMaxHit){

	Distance = glob->TELCorrIntervalY[n];

	for(while2 = 0; while2 < VecDigY[n+1].size(); while2++){

	  if(TMath::Abs(DigY[n][while1] - VecDigY[n+1][while2]) < Distance){
	    Distance = TMath::Abs(DigY[n][while1] - VecDigY[n+1][while2]);
	    DigY[n+1][while1] = VecDigY[n+1][while2];
	    CogY[n+1][while1] = VecCogY[n+1][while2];
	  }

	}//while2
	
	while1++;
	
      }//while1
      
    }//n

    for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){
      
      VecDigX[n].clear();
      VecDigY[n].clear();
      VecCogX[n].clear();
      VecCogY[n].clear();

    }

    for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){  

      while1 = 0;
      while(DigX[n][while1] > 0.0){
	h20[0][n]->Fill(DigX[n][while1],DigX[n+1][while1]);	
	//h20[0][n]->Fill(CogX[n][while1],CogX[n+1][while1]);	
	while1++;
      }

      while1 = 0;
      while(DigY[n][while1] > 0.0){
	h20[1][n]->Fill(DigY[n][while1],DigY[n+1][while1]);	
	//h20[1][n]->Fill(CogY[n][while1],CogY[n+1][while1]);	
	while1++;
      }

    }//n

    TrigTime[0] = Trigger[0];
    TrigTime[1] = Trigger[1];

    Tracks->Fill();

  }//i

  trackfile->cd();
  Tracks->Write();
  trackfile->Close();

  TF1* etaLine = new TF1("etaLine","1.0*x + 0.0",0,1);
  etaLine->SetLineWidth(2);
  etaLine->SetLineColor(2);

  for(Int_t n = 0; n < glob->TELNumberOfPlanes; n++){

    canv10->cd(n+1);
    h10[0][n]->Draw();

    canv10->cd(n + 1  + glob->TELNumberOfPlanes);
    h10[1][n]->Draw();

    canv11->cd(n+1);
    canv11->GetPad(n+1)->SetGridx();
    canv11->GetPad(n+1)->SetGridy();
    h11[0][n]->Draw();
    etaLine->Draw("SAME");

    canv11->cd(n + 1  + glob->TELNumberOfPlanes);
    canv11->GetPad(n + 1 + glob->TELNumberOfPlanes)->SetGridx();
    canv11->GetPad(n + 1 + glob->TELNumberOfPlanes)->SetGridy();
    h11[1][n]->Draw();
    etaLine->Draw("SAME");

  }//n

  /*
  TCanvas* PHDCanv1 = new TCanvas("PHDcanv1","",800,800);
  PHDCanv1->SetFillColor(0);

  PHDCanv1->cd();
  h10[0][0]->Draw();

  TCanvas* PHDCanv2 = new TCanvas("PHDcanv2","",800,800);
  PHDCanv2->SetFillColor(0);

  PHDCanv2->cd();
  PHDCanv2->SetGridx();
  PHDCanv2->SetGridy();
  h11[0][0]->Draw();
  etaLine->Draw("SAME");
  */

  for(Int_t n = 0; n < glob->TELNumberOfPlanes - 1; n++){

    canv20->cd(n+1);
    canv20->GetPad(n+1)->SetGridx();
    canv20->GetPad(n+1)->SetGridy();
    h20[0][n]->Draw("COLZ");
    canv20->cd(n + glob->TELNumberOfPlanes);
    canv20->GetPad(n + glob->TELNumberOfPlanes)->SetGridx();
    canv20->GetPad(n + glob->TELNumberOfPlanes)->SetGridy();
    h20[1][n]->Draw("COLZ");

  }
  
}//HitNtuplWithEta

TELFile::~TELFile(){

  fclose(SiFile);

  delete [] TELRoot;

  delete [] PlaneXADC;
  delete [] PlaneYADC;

  delete [] APara;

}
